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FACTORIZABLE SCHEMES FOR THE EQUATIONS OF FLUID FLOW 


DAVID SIDILKOVER* 


Abstract. We present an upwind high* resolution factorizable (UHF) discrete scheme for the compressible 
Euler equations that allows to distinguish between full- potential and advection factors at the discrete level. 
The scheme approximates equations in their general conservative form and is related to the family of genuinely 
multidimensional upwind schemes developed previously and demonstrated to have good shock-capturing 
capabilities. A unique property of this scheme is that in addition to the aforementioned features it is also 
factorizable , i.e. it allows to distinguish between full- potential and advection factors at the discrete level. 
The latter property facilitates the construction of optimally efficient multigrid solvers. This is done through 
a relaxation procedure that utilizes the factorizability property. 

Key words, incompressible and compressible flow, factorizable schemes, genuinely multidimensional 
upwind schemes, optimal multigrid efficiency 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. The standard numerical methods used for incompressible flow computations do not 
have much in common with the standard methods for compressible flow. The explanation of this fact is 
because difficulties encountered initially in each case were of a very different nature. We shall begin with 
briefly describing these difficulties and the evolution of the numerical methods for the two classes of problems. 

1.1. Dimension- by-dimension methods (compressible flow). The main difficulty encountered 
when constructing numerical methods for compressible flow equations was the possible presence of disconti- 
nuities in the solutions. It took a prolonged effort of numerous researchers to devise what we call now the 
shock-capturing methodology. One of the basic ingredients of the shock- capturing schemes is the so-called 
(approximate) Riemann solvers (or their alternative - flux-splitting techniques) which are used to construct 
a first order accurate scheme. Another important ingredient is the so-called high- resolution mechanism, that 
allows one to combine higher order accuracy with shock- capturing capabilities, i.e., to circumvent Godunov’s 
theorem. It appears in the form of interpolation (or extrapolation) based on a certain smoothness monitor, 
that is usually implemented in the form of a flux- limiter. Many of the relevant issues could be studied on 
one-dimensional (unsteady) model problems. 

The methods developed for one-dimensional problems were extended later to multidimensions in the 
most straightforward way on a dimension- by- dimension basis. 

One practical need was to compute steady flow, both external and internal, in multidimensions. Steady 
problems are generally solved through (pseudo-) time evolution. In other words, the problems are treated in 
a hyperbolic (with respect to time) sense. Multigrid methods became widely used as a mean to accelerate 
convergence to the steady-state. The key ingredient of a multigrid algorithm is the smoother, i.e., a relaxation 
procedure that efficiently reduces the high-frequency error content. The difficulty, however, is that the 
high-frequency error content may be nearly invisible in the residuals (poor measure of h-ellipticity ) of high- 
resolution schemes constructed on a dimension- by- dimension basis. This makes it inherently impossible to 
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132C, NASA Langley Research Center, Hampton, VA 23681-2199, e*mail: sidilkovflicase.edu. 
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construct an efficient smoother. 

Another major difficulty is that, in general, no distinction can be made between the different co-factors 
at the discrete level, i.e., the standard high- resolution schemes are not factorizable. This deficiency of the 
standard discretizations also contributes to the poor performance of standard multigrid solvers. It also leads 
to the loss of accuracy and deteriorating computational efficiency in the case of low speed compressible flow. 

Discrete schemes that are only partially based on the dimension- by- dimension approach were presented 
in [6], [7], [9] and [13]. In these schemes, some of the second order corrections appear in the form of terms 
approximating mixed derivatives. The remaining corrections are based on dimension- by- dimension high- 
resolution mechanism. Schemes of this type when tuned for steady-state computations would have a better 
ft-ellipticity property than the standard ones. This direction, however, has not been explored. 

1.2. Projection-type methods (incompressible flow). Incompressible flow in one dimension is 
trivial. Therefore, numerical analysts had to address the multidimensional problems right away. Incom- 
pressible flow problems (both steady and unsteady) contain an elliptic part, which rules out the use of an 
explicit scheme. 

One of the earlier approaches was to substitute the continuity equations by the pressure Poisson equa- 
tion, derived using the momentum and the continuity equations (see a summary in [10]). By using such a 
representation, decoupling of the elliptic factor (in the form of the pressure Poisson equation) from the rest 
of the system can be achieved. This formulation has been used in conjunction with multigrid as well. In [20], 
it was demonstrated that by deriving and treating the boundary conditions for the pressure properly one 
can obtain the ideal multigrid efficiency both for inviscid and viscous cases. Some different versions of the 
algorithm are discussed in [15] and in [14], where some steps in deriving the pressure Poisson equation are 
made at the discrete level. Two variants of the algorithm are being pursued: one for structured body-fitted 
grids, another for unstructured grids (see [14]). Ideal multigrid efficiency was demonstrated for both on 
several test cases. A possible extension of this approach to the compressible subsonic case is discussed in 
[15]. 

Several methods exist based on discretizing the equations of incompressible flow in their usual (primitive) 
form. It is interesting to note that one of the earlier methods (MAC) (see [8]) relies on a staggered grid 
discretization and uses a pressure Poisson equation derived at the discrete level to update the pressure. 

A known property of a vector field is that it can be represented as a sum of its irrotational and solenoidal 
components. The attractive feature of the staggered grid discretization is that this property can be imitated 
on the discrete level (second-order accurate approximation to the Cauchy- Riemann equations) without pro- 
ducing an odd-even instability. It is problematic, though, to achieve the same on non-st aggered grids. 

When solving the equations of incompressible flow in their primitive form, the important part of the 
process is to satisfy the continuity equation, i.e., to project the velocity field onto the divergence- free manifold 
or, in other words, to discard the irrotational component of the velocity field. Such a process is the main part 
of the projection method [5], whose original version was based on non-staggered grids. It was reformulated 
later using staggered grids. A convenient way to perform the projection step is to introduce an auxiliary 
variable potential and to solve the resulting Poisson equation. This was first done in [1]. 

A commonly used algorithm for solving the incompressible flow equations is the SIMPLE algorithm by 
Patankar & Spalding [12]. It is interesting to note that this algorithm was extended to the compressible 
(subsonic) case (see [11], [2]). 

Applying multigrid methods for solving the fluid flow problems is one of the main subjects of the 
landmark work [3]. Although a few various possibilities are mentioned, the approach, systematically studied 
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and advocated, is based on the staggered-grid discretization and an auxiliary potential variable. It was shown 
that the proper treatment (Distributive Gauss-Seidcl relaxation - shortly DGS) results in the decoupling of 
different co-factors of the system. It was demonstrated later in [4] that an optimal MG efficiency can be 
obtained for the incompressible Euler equations. Some ways of generalizing these ideas to the compressible 
flow were sketched out in [3] as well. 

The Vorticity-Potential (or vorticity and streamfunction in two dimensions) formulation of the flow 
equations seems very attractive, since the different co- factors of the equations decouple. The well known 
difficulty associated with this formulation, though, is the derivation and treatment of the boundary conditions 
for vorticity. Another difficulty associated with this approach is that the components of the vorticity vector 
in three dimensions need to satisfy a certain compatibility condition. Also, it is highly problematic to obtain 
a numerical scheme with shock-capturing capabilities using the Vorticity-Potential formulation. 

The ideal multigrid efficiency for the compressible Euler equations in the subsonic case was demonstrated 
first in [23] using the canonical variables formulation of the equations [22]. However, this formulation cannot 
be generalized to viscous and unsteady cases. 

1.3. About this paper. In order to achieve optimal multigrid efficiency, we have to use projection- 
type methods to solve the discrete equations, i.e., methods that distinguish between the different co- factors 
of the equations. 

The attempts ([2], [23]) to apply projection- type methods to solve the compressible flow equations were 
limited to the subsonic case, since the discretizations used have no shock-capturing capabilities. On the other 
hand the projection-type methods cannot be applied in conjunction with the standard upwind dimension- 
by-dimension discretizations, since the latter are not factorizable (see [15]), i.e., no distinction between the 
different factors can be made at the discrete level. 

There is a need for a factorizable discretization scheme with good shock-capturing capabilities. The 
scheme should rely on a high resolution mechanism and it should be sensitive to the high-frequency er- 
ror content (h- elliptic) at the same time. The flow equations should also be approximated in their usual 
conservative form, since using special formulations leads to the loss of generality. The separation between 
the treatment of co-factors should be achieved through the relaxation procedure, which may rely on some 
special auxiliary variables and/or special forms of the equations. It may look as if too many (seemingly 
contradictory) requirements need to be satisfied by a single discretization. The purpose of this paper is to 
describe a construction of such a discrete scheme. 

A search for a genuinely multidimensional upwind scheme was motivated initially by the necessity to 
devise a discretization that has the high- resolution and h-ellipticity property at the same time. Such a scheme 
for scalar advection problems was proposed in [16], [21]. This approach was extended to the Euler equations 
in [18, 17, 19] in the context of unstructured triangular grids. It was demonstrated that a simple pointwise 
Gauss-Seidel relaxation is stable when applied directly to the high-resolution discrete equation (contrary to 
the case of the standard schemes) and has good smoothing properties. The latter is a manifestation of the 
h-ellipticity property of the discretization. 

It was pointed out in [15] that the genuinely multidimensional approach towards discretizing the com- 
pressible flow equations leads to schemes that are, in addition to other desirable properties, factorizable . 
In this paper we describe such a scheme, i.e. a scheme which is Upwind, High-resolution and Factorizable 
(UHF). Also, we develop a relaxation procedure that utilizes the factorizability property of the scheme. 

2. Preliminaries. In this section we briefly review the genuinely two-dimensional advection scheme, 
introduce some basic properties of the Euler equations and discuss some standard discretization schemes for 
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Fig. 1. Computational grid segment and a control volume. 


the fluid flow equations. 

2,1. Scalar advection schemes. Consider a steady-state conservation law in two dimensions 

-eAu + f(u) x + g(u) y = 0 (2.1) 

where e > 0 is an infinitesimally small number. A general conservative discretization of (2.1) is given by 

M(/+ ~ /-) + ( 9 4- ~ 5-)] — 0 (2.2) 

where f±,g± are the numerical fluxes. A particular discretization scheme for (2.1) can be given by defining 
these numerical fluxes. 

A numerical scheme can be written in the form 

u 0 = J2 c iUi. ( 2 . 3 ) 

i 

Definition 2.1. The scheme is said to be of the positive type if Ci > 0. 

We shall illustrate the construction of some discretization scheme on the linear constant coefficient 
version of (2.1) 


—eAu + au x + bu y = 0. 

A first order upwind scheme can be given by the following numerical fluxes 

f- ~ 5 a (^o + n 3 ) - ^|a|(u 0 - u 3 ) 

9- ~ + u$) — ^|6|(u 0 — U 5 ). 

Assume for simplicity that 


(2.4) 


(2.5) 


b > a > 0. 

A second order scheme with a genuinely multidimensional flavor is given by 

/- = /«-I&(ti 3 -ii 4 ) 

9- — 9- ~ \cl(us - u 4 ). 

This scheme, however, is not of the positive type. 


4 



In order to combine the positivity property with the second order accuracy we need to incorporate the 
second order correction in a nonlinear fashion. A “straightforward” way of doing this is to modify the fluxes 
in the following way 


/- = /“- ^b( u 3 - 

9- — 9- - \ a{u 5 - U4)ip(S _), 

which gives a positive second order scheme if the limiter function ip satisfies the inequality 

ib(Q) 

o < ip(Q) < l; 0<™<i, 


(2.7) 


( 2 . 8 ) 


V>(1) = l 


and the arguments of the limiter function are defined by 


aju o - U3) 
b(u 3 — 114 ) 


(2.9) 


_ ^6(uoj^U5) 
a(u 5 - 114 ) 

Using the following identities 

b(u 3 — u 4 )t/>(Q_ ) = a(u 0 — u 3 )ip(Q — )/Q— ^2 

a(u 5 - U 4 )ip(S~) = -b(uo - U 5 )ip(S-)/R-, 

it is easy to see that this scheme is no better than the standard high-resolution schemes, since under certain 
circumstances it becomes identical to a central scheme (no h-ellipticity property). The fact that it is positive 
for a non-compressive limiter only (i.e., no artificial compression can be added) is a disadvantage as well. 

Note, that a part of the second order correction can be added with no limiter while still resulting in a 
positive (first order) scheme - the so-called A r -scheme (see [16, 21]) 


/- = f- ~ 5 a ( u 3 - «4) (2.11) 

9 - =9-~ | a ( w 5 - u 4 ) 

A second order high-resolution scheme can be obtained by adding the remaining part of the correction in a 
nonlinear fashion 


/_ = f- ~ %(b - a)(u 3 - u 4 )ip(ft-) 


(2.12) 


where 


q(u 5 - u 4 ) 
b(u 3 - 11 4 ) 


It is easy to see ([21]) that (2.12) defines a positive scheme if 


0 < ip{R) < 2; 


0 < 


iP(R) 


- R 


< 2 . 


(2.13) 
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2.2. Euler equations. The quasilinear form of the steady Euler equations in two dimensions can be 
written in the following way 


Lu = 0 (2.14) 

where u (s, u, u, p) is the vector of unknown functions (namely, entropy, horizontal velocity component, 
vertical velocity component and pressure, respectively), 

( q 0 0 

L = 0 pq o 

0 0 pq 

\ 0 pc 2 d x pc 2 d y 

p stands for density, c for the speed of sound and q denotes the advection operator (with velocity (u, u)). We 
shall consider the case of subsonic flow (u 2 + v 2 < c 2 ) throughout the paper. It is sufficient for the purpose 
of the analysis in this paper to consider the linear constant coefficient case. 

Determining the type of a system of partial differential equations can be done formally by computing 
the determinant of matrix L and examining its principal part 


° ) 

0 . 


det L= {pq) 2 ■ ( q 2 - c 2 A). (2.16) 

The first multiplier in (2.16) is the advection operator (times the density) in the power equal to the dimension 
of the problem. The second multiplier is the full-potential operator, which, in subsonic case, is of elliptic 
type. Thus, the steady Euler system of equations is of the mixed hyperbolic-elliptic type. 

2.2.1. One-dimensional case. Consider a first order upwind scheme for the one-dimensional Euler 
equations. Without loss of generality we consider the primitive variable formulation 


Lu = 0 

where u = ( s,u,p) T and 

( ud x 0 0 \ 

0 pud x d x I 
0 pc 2 d x ud x j 

A first order upwind scheme approximating (2.17) is given by 


(2.17) 


(2.18) 


L h u h = 0 (2.19) 

where the discrete variables 

u h = {s h ,u h ,p h ) T (2.20) 

and 

/ -fM^x + Q 2 ' 1 o o 

Lh =\ 0 p(- \cd h xx + Q 2h ) df - 

V 0 pc 2 {d 2h - -§ cdxx + Q2» 

where h is a mesh size, d£ x is a central approximation of the second derivative, d 2h is a central approximation 
of the first derivative and Q 2h = ud 2h is the advection operator. 
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Factorization. The determinant of L h : 


det (L h ) = p{-\\u \d h xx + <? 2h )[(c 2 - u*)d h xx ] (2.22) 

The first factor is the upwind scheme approximating the advection operator corresponding to the entropy 
equation. The Full-Potential factor is approximated by a “short” central difference. The issue of factorization 
appears to be trivial in this case, since the momentum and the pressure equations correspond solely to the 
elliptic factor. 

2.2.2. Two-dimensional case. Since the equation for entropy decouples in the form of the advection 
equation, whose numerical treatment is straightforward, we can consider without loss of generality the case of 
isentropic flow. We can also assume for simplicity that p - 1. Two-dimensional steady isentropic linearized 
Euler equations are now given by the following 

Lu = 0, (2*23) 

where u = (u,u,p) T and 

/ q 0 d x \ 

L = j 0 q d y . (2.24) 

\ (?d x c 2 d y q / 


The FDA (First Differential Approximation or Modified equations) corresponding to the first order upwind 
scheme, constructed on the dimension- by-dimension basis, is given by 


L[ da 


( (q — |(c - |u|)d IX ) 0 d x — 2c^ xx \ 

0 (q — \ ( c ~ M)^yy) ®y~2c^yy I 

V ^(*.-3?®-.) Q-IcA ) 


(2.25) 


where q is an FDA of the first order upwind advection scheme. It is easy to verify that L[ da does not 
preserve the factorizability property of L, i.e., in the det(Lj ?jD ^) we cannot distinguish between the FDAs 
corresponding to advection and full-potential factors. The same applies, of course, to the corresponding 
discrete scheme. 

The dimension- by-dimension scheme (whose FDA is given by (2.25)) can be upgraded to higher order 
accuracy by introducing the standard nonlinear high-resolution correction. However, by doing so, we obtain a 
scheme that is not only not factorizable but also can be insensitive to some high-frequency error components 
(has a poor measure of h-ellipticity). All this makes it practically impossible to construct an efficient solver 
for the resulting discrete equations. 


3. Factorizable scheme. In order to combine the high- resolution and /i-ellipticity properties the ex- 
tension of the ideas leading to the genuinely multidimensional advection scheme need to be generalized to 
systems of equations. Such a generalization on triangular grids was proposed in [18, 17, 19]. As well as 
in the scalar case, the second-order corrections are given by the terms approximating mixed derivatives. 
These corrections may rely on two-dimensional limiters. It was demonstrated that the Gauss-Seidel relax- 
ation is stable when applied directly to the high-resolution discrete equations (unlike in the case of standard 
dimension- by-dimension schemes). This indicated that these new schemes have stability properties superior 
to the standard ones. 

It was pointed out in [17] that some of these corrections can be added without limiters, in a linear 
fashion. Indeed, similar to the scalar case (see [21] and also §2.1), in order to obtain a scheme with better 
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linear stability properties, it is desirable to include as much of the second order correction as possible without 
limiters. This is especially important in the context of structured grids (“9-point box” stencil); otherwise 
one can obtain a scheme that is no better than the standard ones (see §2.1). 

However, unlike the scalar case, it is unclear what is a suitable criterion for determining which parts 
of the corrections do not have to rely on the limiters. The difficulty here is that the notions of maximum 
principle and positivity do not extend to systems of equations in a suitable way. 

We recall here the observation that was made in [15]. Consider the FDA corresponding to the second 
order genuinely multidimensional upwind scheme with all the second order corrections added in a linear 
fashion 

*1 — 2 ” \ u \)&xx ~^( c ~~ \ v \)&xy &x ~ ^cQx \ 

-%{c-\u\)d xy q — §( c ~ \ v \)dyy dy-\\Qv • (3.1) 

<?d x - ±cQ x c 2 d y - | cQy Q - |cA / 

It can be verified that (3.1) is factorizable. 

The implication of this is that the multidimensional corrections not only lead to a second-order scheme, 
but some of them are also responsible for the factorizability. Since it is a desirable property, the latter 
corrections should be included without limiters. 

There exist many discretizations that correspond to (3.1). Note that not all of them are factorizable. 
We would now like to construct an actual discrete scheme which is factorizable as well. 

The technical difficulty here is that the following differential equalities 


rFDA 


&xx&yy &xy&xy 

&xx&y — & X y&x (3.2) 

&yy&x — &xy&y 

need to have discrete analogs. Approximating the derivatives d x ,d xx ,d y ,d yy by standard central finite 
differences and discretizing the mixed derivatives in some reasonable way, we can see that the discrete 
analog of (3.2) does not hold. 

Let us introduce some non-standard finite differences 


d x 


1 1 

f 

-1 

0 

1 \ 

— 11 


1 

2 


8 h 


-2 

0 


' v ~ 8h 


0 

0 

0 



-1 

0 

1 1 



-1 

-2 

-i ) 


£U = -- 


1 1 

( 

1 

-2 

1 

) 

= 1 1 

1 

2 


4 


2 

-4 

2 


; vv ~ 4 ^ 

-2 

-4 

' 2 



1 

-2 

1 

> 


l 1 

2 

i / 


Let us denote also 


o _ 1 J_ 
° xy ~ 4 h* 


-1 0 1 

0 0 0 

1 0 -1 


/ 1 2 1 \ 
A = \h 2 - 12 2 

V 1 2 1/ 


(3.3) 


(3-4) 


(3.5) 


(3.6) 
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and 


Qx — “ 1 “ V&xy ^3 7 ) 

Q x = ud X y +vdyy. 

It is easy to see that the discrete analog of (3.2) holds if the derivatives are approximated by the finite 
differences as defined by (3.3), (3.4) and (3.5). 

Consider now the following discrete scheme that corresponds to the FDA given by (3.1) 

/ q h - |(c- \u\)d h xx -*(c- M)3 X \ % - \ 

L h =\ -l(c-\u\)d* y qfc-f (c-\v\)B^ y d h y -±\Q h y • (3-8) 

V c 2 d h x - f cQ* c 2 ^ - | cQ y Q h - \ch h ) 

It is factorizable and its determinant is given by 

det(L h ) = q W - ^A^q* - ^(c- \u\)d* x - ^(c- \v\)B£ y ) 

-(B h x - - jcQ.) 

Note that the discrete approximation of the elliptic factor becomes the so-called “skewed 5 discrete Laplacian 
(not an h-e Uiptic operator) when the flow speed (Mach number) goes to zero. This difficulty can be dealt 
with in various ways. The simplest option is to slightly modify the scheme so that the problem mentioned 
above disappears. 

Denote the total velocity 

U = |\AtM^|, 


the Mach number 



and the “augmented 51 Mach number 


M = 


M, if M > Ch 
Ch otherwise 


(3.9) 


where C > 0 is a constant. 

We can “rescale 55 the artificial dissipation terms, obtaining the following discrete scheme, which we shall 
call UHF (Upwind, High-resolution, Factorizable) 


L = 


( q* - - |«|)a£ x 

-h(U -\u\)d h xy 


q/v _ | ([7 _ \ v]) BH y 

c 2 (dy-I^Qy) 


o/i h M 

U x 2 c 


Qt 

h m s\h 


fsh h M 

°y 2 c^V 
h _ h c Ah 
2 M 


Q 


(3.10) 


where 


= q h - -h 2 cCA b , 
2 


(3.11) 
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is the discrete advection operator augmented by a second-order small dissipative term and A h is the usual 
“5-point-star” discrete Laplacian. This “regularization,” obviously, prevents the advection scheme from 
degenerating at a stagnation point. It plays an additional important role, though, which will be explained 
later in this section. 

Before discussing some properties of the UHF scheme and computing its determinant, we shall introduce 
the simplified notation 



Q h =d h - ^~^-O h 

* * 2 cM 

d h = d h - ~^~o h 

v u y 2 cM 


(3.12) 


Q h = Q h - 

X u x 2 c 

d h = d h - 

V V 2c 


(3.13) 

and 



(3.14) 


q 1 ' = q k - - M01 + (u- 


(3.15) 

Then the UHF scheme (3.10) 

can be written as 





9* \ 


L = 

-%{U-\u\)d h xy q*-| (U-)v\)d£ y 

K 

(3.16) 


c 2 d h x <?8 h y 

q* J 


and its determinant 

det {L h ) = q* ■ [c 2 (d h x d h x + d h y 8 h y ) - q h • q'*]. 


(3.17) 

Let us take a closer look at what happens at a stagnation point (U = u = v 

= 0) 



q h = —-h 2 cCA h 
2 


(3.18) 


q* = ~h 2 cCA h 


(3.19) 


q h = — —A h . 


(3.20) 

It becomes plain that 

q h ■ q h = i/i 2 c 2 A h A h , 


(3.21) 


which is an /i-elliptic operator. We can conclude now that for the UHF scheme (3.16) (or (3.10)) the 
approximation of the Full- Potential factor remains /i-elliptic in the case of vanishing velocity. Moreover, 
it can be verified that the Full- Potential factor is approximated in this case by a “regular” five-point star 
discrete Laplacian A*. 
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4. The relaxation procedure. Having constructed a factorizable discretization scheme for the Euler 
equations, an important question becomes how to use this property in order to obtain an efficient algorithm 
for solving the discrete equations. 

Introduce the auxiliary variables, namely, the discrete stream-function and potential by the following 

/ u h 

1 / 1IFM* 1 

(4.1) 




where 


v P h 


M h = 




( d y B h x \ 

-d h x d h y 

V 0 q h / 


We can also define the discrete vorticity 


= curl ' • (u h ,v h ) T = (& n y , -8' x ) ■ (u h , v h ) T 


(4.2) 


(4.3) 


or 


w h = (d h x d h x + d h y d h y )* h . 


(4.4) 


Introduce also the following matrix operator 



(4.5) 


It is easy to see that applying the discrete curl /> operator to the momentum equations (pre-multiplying L h 
by V h ) and performing the substitution of variables according to (4.1) (post-multiplying L h by M h ), wc 
obtain 


V h ■ L h ■ M h 


q h ■ (d*d x + o*dy) 


ha fl \ 


<?(d h x d h x + d*d n y )~ q h q 


A 


~h -h 


(4.6) 


In other words we end up with “solving 15 the system 



The Distributive relaxation procedure at a point amounts to computing updates to the discrete potential 
and vorticity (or streamfunction) according to (4.8) and to translating them into corresponding updates of 
the velocity components and the pressure at the point of interest and its neighbors according to matrix (4.2) 
and (4.1). 

Note that the operator {d h x d h x + &ydy) is not ^-elliptic. Therefore, the relaxation procedure described 
above will not smooth certain high-frequency error components. This lack of /i-ellipticity seems unavoidable 
for obtaining the desired factorization. However, a simple remedy exits for this trouble. It is to augment 
each sweep of the Distributive relaxation by a sweep of a point Collective Gauss-Seidel relaxation. The latter 
will smooth the problematic error components. 

It might also be useful to perform the point collective relaxation, instead of Distributive, on and near 
the boundary, thus avoiding the difficulty of imposing boundary conditions for vorticity. 
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5. Discussion and conclusions. When discretizing a system of partial differential equations, being 
concerned just with obtaining an approximation of a certain order of accuracy may not necessarily lead 
to a satisfactory result. It may be very useful to make sure that the discretization also imitates some 
fundamental properties of the PDEs. We have outlined a construction of a discretization scheme that not 
only approximates the compressible Euler equations, but imitates their factor izability property. This paves 
the way towards construction of optimally efficient multigrid solvers and also should alleviate the problem 
associated with computation of low-speed flow using the standard shock- capturing schemes. 

This paper addressed the subsonic case only. The constructed UHF scheme, though, is related to the 
family of the genuinely multidimensional upwind schemes constructed previously (see [18, 17, 19]). Schemes 
of this type were demonstrated to have excellent shock capturing capabilities in transonic-supersonic cases 
together with maintaining the /i-ellipticity property. Therefore, we do not anticipate any difficulties in 
extending the factorizable scheme to the transonic- supersonic regimes. This is a subject of future work. 

We also developed a relaxation procedure that uses the factorizability property of the scheme for the 
purpose of obtaining ideal multigrid efficiency. This procedure relies on auxiliary potential and vorticity 
(or streamfunction) variables. Note that using only an auxiliary potential variable (as in the incompressible 
case) is not sufficient yet to decouple the different factors (or obtaining upper- or lower-triangular matrix 
of difference operators). This is due to the bulk viscosity- like terms, which are a part of the artificial 
dissipation in the constructed scheme and are essential for its factorizability and second order accuracy. 
Another auxiliary variable is needed - vorticity (or streamfunction). 

The approach we introduce in this paper seems very general. The future work will be devoted generalizing 
it to three dimensions, to Navier-Stokes equations and to unsteady problems. 

The proposed approach also essentially unifies the numerical treatment of the steady incompressible and 
compressible flow, since, on one hand, it belongs to the class of projection- type methods, on the other hand 
its extension to the supersonic case has shock- capturing capabilities. 
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